rm(list=ls())

library(rstan) 
library(ggplot2)
library(Hmisc)
library(foreign)
library(car)
library(lme4) 
library(shinystan)
library(readr)

multi_final<-  read_csv("(...)/pseudo_final.csv")
attach(multi_final)


# Create a non-missing dataset for education

subset_edu <- c("cnt", "av_spend_edu","std_av_party", "std_av_income", "std_gini_net", "std_unemp", "std_gdp_growth","std_lagedu_av","groupid", "year")
multi_final_edu <- multi_final[subset_edu]     
#View(multi_final_edu)


multi_final_edu2 <-   na.omit(multi_final_edu)
#View(multi_final_edu2)
attach(multi_final_edu2)

describe(cnt)
cnt <- recode(cnt,"1=1;2=2;3=3;4=4;5=5;6=6;7=7;8=8;9=9;10=10;11=11;12=12;13=13;14=14;15=15;16=16;17=17;18=18;20=19;21=20")

#Prepare the variables for the analysis

recyear <- as.numeric(year == 2006)
y <- multi_final_edu2$av_spend_edu
summary(y)
groupid1 <- as.numeric(groupid == 1)
groupid2 <- as.numeric(groupid == 2)
groupid3 <- as.numeric(groupid == 3)
groupid4 <- as.numeric(groupid == 4)
groupid5 <- as.numeric(groupid == 5)
groupid6 <- as.numeric(groupid == 6)
groupid7 <- as.numeric(groupid == 7)
groupid8 <- as.numeric(groupid == 8)
groupid9 <- as.numeric(groupid == 9)
groupid10 <- as.numeric(groupid == 10)
groupid11 <- as.numeric(groupid == 11)
groupid12 <- as.numeric(groupid == 12)
groupid13 <- as.numeric(groupid == 13)
groupid14 <- as.numeric(groupid == 14)
groupid15 <- as.numeric(groupid == 15)
groupid16 <- as.numeric(groupid == 16)
groupid17 <- as.numeric(groupid == 17)
groupid18 <- as.numeric(groupid == 18)
groupid19 <- as.numeric(groupid == 19)
groupid20 <- as.numeric(groupid == 20)
groupid21 <- as.numeric(groupid == 21)
groupid22 <- as.numeric(groupid == 22)
groupid23 <- as.numeric(groupid == 23)
groupid24 <- as.numeric(groupid == 24)
groupid25 <- as.numeric(groupid == 25)
groupid26 <- as.numeric(groupid == 26)
groupid27 <- as.numeric(groupid == 27)

party <- round(std_av_party,digits=2) 
income <- round(std_av_income,digits=2)
gini <- round(std_gini_net,digits=2)
unemployment <- round(std_unemp,digits=2)
lagedu_av <- round(std_lagedu_av,digits=2)
growth <- round(std_gdp_growth,digits=2)
J<-20
describe(cnt)

attach(multi_final_edu2)


# 1st part of the Bayesian model, define the data
modeledu_av = "
data {
int<lower=0> N; 
int<lower=1> J; 
int<lower=1,upper=J> cnt[N];
real <lower=-4.35,upper=4.44> party[N];
real <lower= -2.46,upper=2.81 >income[N];
real <lower=-1.59,upper=2.02> gini[N];
real <lower=-1.29,upper=2.52 > unemployment[N];
real<lower= -1.87,upper=3.04 > growth[N];
real<lower=  -2.03,upper=2.19> lagedu_av[N];
int<lower=0,upper=1> recyear[N];
int<lower=0,upper=1> groupid2[N];
int<lower=0,upper=1> groupid3[N];
int<lower=0,upper=1> groupid4[N];
int<lower=0,upper=1> groupid5[N];
int<lower=0,upper=1> groupid6[N];
int<lower=0,upper=1> groupid7[N];
int<lower=0,upper=1> groupid8[N];
int<lower=0,upper=1> groupid9[N];
int<lower=0,upper=1> groupid10[N];
int<lower=0,upper=1> groupid11[N];
int<lower=0,upper=1> groupid12[N];
int<lower=0,upper=1> groupid13[N];
int<lower=0,upper=1> groupid14[N];
int<lower=0,upper=1> groupid15[N];
int<lower=0,upper=1> groupid16[N];
int<lower=0,upper=1> groupid17[N];
int<lower=0,upper=1> groupid18[N];
int<lower=0,upper=1> groupid19[N];
int<lower=0,upper=1> groupid20[N];
int<lower=0,upper=1> groupid21[N];
int<lower=0,upper=1> groupid22[N];
int<lower=0,upper=1> groupid23[N];
int<lower=0,upper=1> groupid24[N];
int<lower=0,upper=1> groupid25[N];
int<lower=0,upper=1> groupid26[N];
int<lower=0,upper=1> groupid27[N];
vector[N] y;
} 

parameters {
vector[J] a;
vector[33] beta;
real mu_a;
real<lower=0> sigma_a;
real<lower=0> sigma_y;
}



transformed parameters {

vector[N] y_hat;
for (i in 1:N)
y_hat[i] <- a[cnt[i]] +beta[1]*party[i]+ beta[2]*income[i] + beta[3]*gini[i]+ beta[4]*unemployment[i] + beta[5]*growth[i]+beta[6]*lagedu_av[i]+
beta[7]*recyear[i]+beta[8]*groupid2[i]+beta[9]*groupid3[i]+
beta[10]*groupid4[i]+beta[11]*groupid5[i]+beta[12]*groupid6[i]+beta[13]*groupid7[i]+beta[14]*groupid8[i]+beta[15]*groupid9[i]+
beta[16]*groupid10[i]+beta[17]*groupid11[i]+beta[18]*groupid12[i]+
beta[19]*groupid13[i]+beta[20]*groupid14[i]+beta[21]*groupid15[i]+beta[22]*groupid16[i]+beta[23]*groupid17[i]+beta[24]*groupid18[i]+beta[25]*groupid19[i]+beta[26]*groupid20[i]+
beta[27]*groupid21[i]+beta[28]*groupid22[i]+beta[29]*groupid23[i]+beta[30]*groupid24[i]+beta[31]*groupid25[i]+beta[32]*groupid26[i]+beta[33]*groupid27[i];
}


model {
mu_a ~ normal(0, 10);
a ~ normal(mu_a, sigma_a);
beta ~ normal(0, 10);
y ~ normal(y_hat, sigma_y);
sigma_a ~ cauchy(0, 2.5);
sigma_y ~ cauchy(0, 2.5);
}

"
sm_edu_av <- stan_model(model_code= modeledu_av)


dataList.1 <- list(N=length(y),y=y,J=J, cnt=cnt,party=party,income=income,gini=gini, unemployment=unemployment,growth=growth,lagedu_av=lagedu_av, recyear=recyear,
                   groupid1=groupid1,groupid2=groupid2,groupid3=groupid3,groupid4=groupid4,groupid5=groupid5,groupid6=groupid6,groupid7=groupid7,groupid8=groupid8,groupid9=groupid9,
                   groupid10=groupid10,groupid11=groupid11,groupid12=groupid12,groupid13=groupid13,groupid14=groupid14,groupid15=groupid15,groupid16=groupid16,groupid17=groupid17,groupid18=groupid18,groupid19=groupid19,
                   groupid20=groupid20,groupid21=groupid21,groupid22=groupid22,groupid23=groupid23,groupid24=groupid24,groupid25=groupid25,groupid26=groupid26,groupid27=groupid27)

fit_edu_av <- sampling(sm_edu_av, data=dataList.1, iter=5000,chains=5)


shinystan_edu_av <- as.shinystan(fit_edu_av)
launch_shinystan(shinystan_edu_av)


